
SUBROUTINE force_erp (mjd,rsat,vsat,r_sun,fx_erp,fy_erp,fz_erp)

! ----------------------------------------------------------------------
! SUBROUTINE: force_erp.f90
! ----------------------------------------------------------------------
! Purpose:
! Acceleration due to the earth radiation pressure. 
! This model is only for GNSS satellites, not for LEO satellites which needs
! more precise model, such as a numerical one.
! ----------------------------------------------------------------------
! Input arguments:
! - mjd          : time variable
! - rsat         : satellite position vector (m)
! - vsat         : satellite velocity vector
! - r_sun        : Sun position vector
! 
! Output arguments:
! - fx,fy,fz:	 : Accelerations in the inertial frame (m/s^2)
! ----------------------------------------------------------------------
! Author :	Dr. Tzupang Tseng
!
! Created:	15-03-2018
!
! Changes:      03-10-2018 Tzupang Tseng: Add the conversion between PRN and SVN
!                                         ADD the antenna thrust info for GNSS 
!
! Copyright:  GEOSCIENCE AUSTRALIA
! ----------------------------------------------------------------------

      USE mdl_precision
      USE mdl_num
      USE mdl_param
      IMPLICIT NONE

! ----------------------------------------------------------------------
! Dummy arguments declaration
! ----------------------------------------------------------------------
      REAL (KIND = prec_q)               :: mjd  
      REAL (KIND = prec_q), DIMENSION(3) :: rsat,vsat,r_sun
      REAL (KIND = prec_q)               :: fx_erp,fy_erp,fz_erp

!! repro3 albedo model
       INTEGER*4 ERM,ANT,GRD,REFF,MONTH,SVN
       REAL*8 KAPPA(3,3)
       REAL (KIND = prec_q), DIMENSION(3) :: SUN,ACCEL
       REAL (KIND = prec_q), DIMENSION(6) :: YSAT

!! Force generated by repro3 routines
! -----------------------------------
ERM = 1
ANT = 0
GRD = 1
REFF = 0 
YSAT(1:3) = rsat
YSAT(4:6) = vsat
KAPPA = 0.d0
MONTH = 1

!print*,'Block type, SVNID, BLKID,  power, MASS = ', BLKTYP, SVNID, BLKID,power, MASS
CALL ERPFBOXW(ERM,ANT,GRD,REFF,YSAT,R_SUN,KAPPA,MONTH,SVNID,MJD,ACCEL)

fx_erp = ACCEL(1)
fy_erp = ACCEL(2)
fz_erp = ACCEL(3)
!IF (BLKTYP(1:3) == 'GAL' .OR. BLKTYP(1:3) == 'GLO') THEN
!PRINT*,'fx_erp, fy_erp, fz_erp =', fx_erp, fy_erp, fz_erp
!END IF

END SUBROUTINE
